Freak Waves in Random Oceanic Sea States 
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Freak waves are very large, rare events in a random ocean wave train. Here we study the 
numerical generation of freak waves in a random sea state characterized by the JONSWAP power 
spectrum. We assume, to cubic order in nonlinearity, that the wave dynamics are governed by the 
nonlinear Schroedinger (NLS) equation. We identify two parameters in the power spectrum that 
control the nonlinear dynamics: the Phillips parameter a and the enhancement coefficient 7. We 
discuss how freak waves in a random sea state are more likely to occur for large values of a and 
7. Our results are supported by extensive numerical simulations of the NLS equation with random 
initial conditions. Comparison with linear simulations are also reported. 
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Freak waves are extraordinarily large water waves whose heights exceed by a factor of 2.2 the significant wave height 
O |. of a measured wave train Q . The mechanism of freak wave generation has become an issue of principal interest due 
' to their potentially devastating effects on offshore structures and ships. In addition to the formation of such waves 
in the presence of strong currents Q or as a result of a simple chance superposition of Fourier modes with coherent 
phases, it has recently been established that the nonlinear Schroedinger (NLS) equation can describe many of the 
features of the dynamics of freak waves which are found to arise as a result of the nonlinear self-focusing phenomena 
[^-^|. The self- focusing effect arises from the Benjamin- Feir instability |j: a monochromatic wave of amplitude, ao, 
f"*^ [ and wave number, ko, modulationally perturbed on a wavelength L = 2tt/ Afc, is unstable whenever Ak/(ko>e) < 2y/2, 
where e is the steepness of the carrier wave defined as £=fco&o- The instability causes a local exponential growth 
in the amplitude of the wave train. This result is established from a linear stability analysis of the NLS equation 
H and has been confirmed, for small values of the steepness, by numerical simulations of the fully nonlinear water 
wave equations (||| (f° r high values of steepness wave breaking, which is clearly not included in the NLS model, 
can occur). Moreover, it is known that small-amplitude instabilities are but a particular case of the much more 
complicated and general analytical solutions of the NLS equation obtained by exploiting its integrability properties 
via Inverse Scattering theory in the ^-function representation |l~HJ[2| ] . 

Even though the above results are well understood and robust from a physical and mathematical point of view, it 
is still unclear how freak waves are generated via the Benjamin- Feir instability in more realistic oceanic conditions, 
i.e. in those characterized not by a simple monochromatic wave perturbed by two small side-bands, but instead 
' by a complex spectrum whose perturbation of the carrier wave cannot be viewed as being small. Furthermore, the 
focus herein is not to attempt to model ocean waves but instead to study leading order effects using the nonlinear 
Schroedinger equation, as suggested by Research at higher order suggests that the results given herein are 

indicative of many physical phenomena in the primitive equations [^|^] . 

In this Letter our attention is focused on freak wave generation in numerical simulations of the NLS equation where 
we assume initial conditions typical of oceanic sea states described by the JONSWAP power spectrum (see, e.g. p3|): 
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where er =0.07 if / < /o and cr =0.09 if / > Jo- Our use of the JONSWAP formula is based upon the established result 
that developing storm dynamics are governed by this spectrum for a range of the parameters [[l3|. The constants a, 7 
and <tq were originally obtained by fitting experimental data from the international JONSWAP experiment conducted 
during 1968-69 in the North Sea. Here /o is the dominant frequency, 7 is the "enhancement" coefficient and a is the 
Phillips parameter. For 7=1 and a=0.008I the spectrum reduces identically to that of Pierson and Moskowitz fTs] ] 
which describes a fully developed sea state, i.e. one which has evolved over infinite time and space. For / = /o, the 
expression E — 7 cx p[ -1 is effectively equal to the value of 7 and, as / moves away from /o in either direction, E tends 
rapidly to unity. Therefore as 7 increases, the spectrum becomes higher and narrower around the spectral peak. In 
Fig. [jjwe show the JONSWAP spectrum for different values of 7 (7=1, 5, 10) for /o=0.1 Hz and a=0.0081. 

The major finding we would like to discuss herein is that as 7 and a grow, the nonlinearity becomes more important 
and the probability of the formation of freak waves increases. Our results have been achieved by considering the NLS 
equation as the simplest nonlinear evolution equation for describing deep-water wave trains. We have performed 
numerical simulations using the JONSWAP spectrum to determine the initial conditions. Since the analytical form of 
the spectrum is given as a function of frequency, the analysis is carried out by considering the so called time- like NLS 
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equation (TNLS) (for the use of time-like equations in water waves see, e.g., [|l4|-|l6|) which describes the evolution of 
the complex envelope A in deep water waves: 

A x + i(^\ A tt +ie 2 | A\ 2 A = 0, (2) 

where dimensional quantities denoted with primes have been scaled according to: A — ciqA' , x — x'/ko and t — t'/Au) 
with 1 / Aui a characteristic time scale of the envelope which corresponds to the width of the frequency spectrum. Eq. 
(||) solves a boundary value problem: given the temporal evolution A(0,t) at some location x — 0, eq. (2) determines 
the wave motion over all space, A(x,t). 

At this point, it is instructive to introduce a parameter that estimates the influence of the nonlinearity in deep 
water waves. This parameter, which is a kind of "Ursell" number can be obtained as the ratio of the nonlinear 
and dispersive terms in the TNLS equation: 

U (aw/w ) ^ 

When Ur << 1 waves are essentially linear and their dynamics can be expressed as a simple superposition of sinusoidal 
waves. For Ur > 1, the dynamics become nonlinear and the evolution of the wave train is likely dominated by envelope 
solitons or unstable mode solutions such as those studied by Yuen and coworkers || . 

Many aspects of the importance of the nonlinearity can be addressed by computing Ur from the spectrum (Q). In 
Fig. [2] we show the Ursell number as a function of the parameter 7 for a=0.0081 and a=0.0162. In the construction 
of the plot an estimation of e and Alu/uio needs to be given. The steepness e has been estimated as the product of 
the wave number, ko, of the carrier wave with a characteristic wave amplitude which we compute as the significant 
wave height (the mean of the highest 1/3 wave heights in a wave train), H s , divided by 2. Am is a measure of the 
width of the spectrum and it has been estimated as the half-width at half-maximum. From the plot it is evident 
that for the Pierson-Moskowitz spectrum (7 = 1) the Ursell number is quite small: this indicates that dispersion 
dominates nonlinearity. Formally, the NLS equation is derived assuming that the spectrum is narrow banded and 
the steepness is small. It has to be pointed out that for small values of 7 (7 = 1,2) the spectrum is not narrow 
banded; as 7 increases the spectrum becomes narrower (Auj/loq ~0.2 or less), suggesting that the NLS equation is 
more appropriate. For large values of 7 the mean steepness increases; for 7 = 8, a = 0.01 the steepness is equal to 
0.16, therefore the equation is no longer valid and higher order terms in steepness are required. In Fig. ^| we have 
placed vertical lines at 7 — 2.5 and 7 = 8 to indicate the region in which the NLS equation is applicable. 

When the spectral width becomes large, one expects results which are somewhat out of the range of applicability of 
the NLS equation. As pointed out in a number of papers p|JTo|] the main defect in the NLS equation concerning the 
narrow-band approximation, arises from the fact that linear dispersion is not at high enough order. Reference |ic| l 
proposes an equation that includes all the terms in the linear dispersion relation. The equation (eq. (1) in their paper), 
which is basically the NLS equation with the full linear dispersion relation of the primitive equations, "reproduces 
exactly the conditions for nonlinear four-wave resonance even for bandwidth greater than unity" . In the linear limit 
the equation is exact. In our numerical simulations with the Pierson Moskowitz spectrum (7 = 1, a = 0.0081 ) we 
have used both the NLS equation and its modified form (eq. (1) in [|l0|]). In this specific case the two equations give 
basically the same results: nonlinearities are weak (Ursell number=0.03, see Fig. |) and the dynamics are basically 
linear: the correction in the linear dispersion relation does not essentially alter the value of the maximum simulated 
wave amplitudes. The results of these tests have convinced us that, for the important range 7 = 2.5 — 8, the simpler 
NLS eqution is a valid approach for studying many of the properties of rogue waves. 

The influence of the parameter a consists in increasing the energy content of the time series and, therefore as a 
increases, the wave amplitude and consequently the wave steepness also increase. If a doubles, the steepness increases 
by a factor of \/2 and the Ursell number by a factor of 2 since the spectral width remains constant. From this analysis 
we expect that large amplitude freak waves (large with respect to their significant height) are more likely to occur 
when 7 and a are both large. 

We now consider numerical simulations of eq. (||) which have been computed using a standard split-step, pseudo- 
spectral Fourier method HI. Initial conditions for the free surface elevation £(0, i) have been constructed as the 
following random process |17|] : 

JV 

C(0,i) = C„ cos(27r/ n t -<M, (4) 



where (f> n are uniformly distributed random numbers on the interval (0,27r), and C n — \j2P(f„)Af n , where P(f) is 
the JONSWAP spectrum given in ([!]). For the computational domain considered, it was checked that the shape of the 
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JONSWAP spectrum has been not altered during the evolution of the NLS equation. It was found that during the 
evolution there is a continuos exchange of energy among the frequencies around the peak. If the numerical time series 
at a fixed spatial location is split into shorter time series in order to compute an average spectrum (this procedure 
is usually adopted when dealing with experimental time series), the shape of the spectrum is preserved with high 
accuracy, at least for 7 < 8. In Fig. || we show an image of smoothed contours of a space-time field of \A\ from a 
numerical simulation of TNLS obtained with 7=4. The dominant frequency and the Phillips parameter of the initial 
wave train were set respectively to 0.1 Hz and a=0.02. A large amplitude wave appears in the simulation and in 
order to better visualize it, in Fig. || we show a time series of the free surface ((t) at x = 1550 m obtained using the 
following relation: 

C{t) = {A(ty 2 * fot + c.c.)/2. (5) 

where c.c. denotes complex conjugate. A freak wave of about 18.5 m in a random wave train with significant wave 
height of H s — 6.9 m is evident at time t = 140 s. We point out that the same simulation, not reported here for 
brevity, with exactly the same initial conditions, has been performed after deleting the nonlinear term in the TNLS 
equation, i.e. the term | A | 2 A was set identically to zero. No waves fulfilling the freak wave threshold (H > 2.2Hs) 
were found. In this linear simulation, the Benjamin- Feir instability cannot occur and "freak" waves can occur only 
via a simple superposition of Fourier modes. This result indicates that, in our simulations, nonlincarity plays an 
important role in the dynamics of freak wave generation. 

In order to give additional quantitative results wc have performed more that 300 simulations of the TNLS equation. 
The simulations have been performed in dimensional units in the following way. An initial time series of 250 seconds 
has been computed from the JONSWAP spectrum for different values of a (from a=0.0081 to a=0.02) and 7 (7=1, 
7=4 and 7=10 ). To increase the number of statistical events we made computer runs with 10 different sets of random 
phases, (f>„. The time series were then evolved according to the TNLS for a distance of 10 km, saving the output 
every 10 m. From an experimental point of view this approach corresponds to setting 1000 probes along the wave 
propagation direction (one every 10 m) and measuring each time series for 250 s at a sampling frequency of 2.05 Hz. 
The significant wave height, H s , of each realization has been computed; the highest wave, H max , has been found and 
the ratio H max /H s has been determined. In order to verify that, in all the simulations performed, our findings are 
really a consequence of the nonlinear dynamics, we have also computed exactly the same simulations using the linear 
version of the TNLS equation. The results are summarized in Fig. ||, |[ Fig. || corresponds to 7 = 1. A horizontal line 
at H max /H s = 2.2 indicates the threshold that arbitrarily discriminates the height of rogue waves. For the Pierson 
Moskowitz spectrum (7 = 1 and a = 0.0081) only one realization of the 10 considered shows a "rogue" wave with 
H m ax/ H s — 2.25. For higher values of a only a few of the realizations show waves with H max /H s slightly greater 
than 2.2. From the plot it is clear that the effects of nonlinearities are rather small for this case (7 = 1). Among 
all the 50 linear simulations performed with 7 = 1, we have encoutered a number of large amplitude waves but none 
exceeds 2.2 H s . 

For 7=4, see Fig. (|6|), the physical picture becomes much more interesting: while in the linear simulations there 
are no freak waves, 50% of the nonlinear simulations performed show at least one freak wave. 

For 7=10 the picture is qualitatively the same and therefore the graphs are not reported. There is clear evidence 
that increasing 7 increases the probability of freak wave occurrences; high values of 7 do not, however, guarantee 
the presence of a giant wave. The local properties of the wave trains are presumably of fundamental importance for 
understanding the formation of freak waves: it may happen that the Benjamin- Feir instability mechanism is satisfied 
only in a small temporal portion of the full wave train, giving rise to a local instability and therefore to the formation 
of a freak wave. 

From a physical point of view, we are aware of the fact that the NLS equation overestimates the region of instability 
and the maximum wave amplitude with respect to higher order models [18|], especially for e greater than 0.1. Further- 
more it is well known that the NLS equation is formally derived from the Euler equations under the assumption of a 
narrow-banded process. Nevertheless, in spite of these deficiencies in the NLS equation, we believe that our results 
provide new important physical insight into the generation of freak waves. Simulations with higher order models |lS| ] 
or directly with the fully nonlinear equations of motion will be required in order to confirm these results. Wave tank 
experiments will also be very useful in this regard. 

Another issue that has to be taken into account for future work is directional spreading: it is well known that sea 
states are not fully unidirectional and directionality can play an important role in the dynamics of ocean waves. In a 
recent paper []l9f we have considered simple initial conditions using the NLS equation in 2+1 dimensions and we have 
found the ubiquitous occurrence of freak waves. Whether the additional directionality in the JONSWAP spectrum 
changes our statistics is still an open question; at the same time we are confident that our results can apply to the 
case in which the spectrum is quasi-unidirectional. In particular, as recently suggested p0| , the so called "energetic 
swells" , which correspond to the early stage of swell development (still characterized by a highly nonlinear regime) , 



3 



are described by high values of 7 and a. Therefore, these sea states are candidates for the occurrence of freak waves. 
In such conditions 1+1 NLS represents a good evolution equation. 
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FIG. 1. The JONSWAP spectrum for 7=1 (dashed line), 7=5 (dotted line), 7=10 (solid line) with /o=0.1 Hz and q=0.0081. 
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FIG. 2. The Ursell number as a function of 7 for the JONSWAP spectrum. 
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FIG. 3. Nonlinear Schroedinger space-time evolution of a random wave field using the JONSWAP spectrum with 7=4, 
a=0.0081. For details refer to the text. Gray scale ranges from m (white) to 12.8 m (dark). 
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FIG. 4. Free surface elevation £(t) at x = 1550 m obtained from Fig. 3. 
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FIG. 5. H ma x/H s as a function of a for 7 = 1. Circles and crosses corresponds respectively to the nonlinear (TNLS) and 
linear simulations. For each value of a, 10 different realizations corresponding to 10 different sets of random numbers have 
been performed. 
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FIG. 6. H max /H s as a function of a for 7 = 4. Circles and crosses corresponds respectively to the nonlinear (TNLS) and 
linear simulations. For each value of a, 10 different realizations corresponding to 10 different sets of random numbers have 
been performed. 
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